% March, 2012
%
% Computes the scattering of a plane wave with one single object
% kh    Helmholtz parameter
% angle Angle of the incoming wave

function oneobj(kh,theta_in)

check_sol = true;
verif2d   = true;

% number of points per wavelength
ppl = 20;
acc = 1e-10;

nobj            = 1;
% The number of points on each contour.
n               = 200;  
ntot            = n*nobj;
% Grid size for 2d visualization
ngrid           = 250;

contours        = cell(nobj,1);
Ctot            = [];
hmax_phy        = 0;

[C,len,xin,tmp] = get_geometry_star(n,[0;0],5,0.);

contours{1}     = C;
h_phy           = max(sqrt(C(2,:).^2 + C(5,:).^2));    
hmax_phy        = max(hmax_phy,h_phy);
Ctot            = [Ctot,C];          

xctot       	= mean(C(1,:));
yctot       	= mean(C(4,:));
Rtot            = max(sqrt((C(1,:) - xctot).^2 + (C(4,:) - yctot).^2 ));

lambda = 2*pi/kh;

if verif2d
  figure(2), hold all
  for i=1:nobj
    plot(contours{i}(1,:),contours{i}(4,:),'k')
  end
  figure(3), hold all
  for i=1:nobj
    plot(contours{i}(1,:),contours{i}(4,:),'k')
  end
end

% Check the size of the discretization vs the wavelength
fprintf(1,'h / lambda = %g\n', kh*hmax_phy/(2*pi))

if kh*hmax_phy/(2*pi) > 1./10
  error('mesh too coarse')
end

if verif2d
  xx      = linspace(xctot - 3*Rtot,xctot + 3*Rtot,ngrid);
  yy      = linspace(yctot - 3*Rtot,yctot + 3*Rtot,ngrid);
  [xxx,yyy] = meshgrid(xx,yy);
  %%% Determine which nodes in the visualization mesh are inside the contour
  %%% by evaluating Cauchy integrals.
  %%% Note that one integral is evaluated for each node in the mesh so this
  %%% can be slightly expensive. 
  %%% If you want a finely resolved visualization mesh, it might be worth
  %%% it to precompute "ind" and store it.
  zzz = [reshape(xxx,1,numel(xxx));...
	 reshape(yyy,1,numel(yyy))];
  nz  = size(zzz,2);
  zz =     (ones(nz,1)*Ctot(1,:) - zzz(1,:)' * ones(1,ntot)) + ...
	   1i * (ones(nz,1)*Ctot(4,:) - zzz(2,:)' * ones(1,ntot));
  h   = 2*pi / n;
  dz  = h*Ctot(2,:)' + h*1i*(Ctot(5,:)');
  uu  = (1/(2*pi*1i))*(1./zz)*dz;
  ind = find(abs(uu) > 0.01);
else
  % Verify the solution on a circle
  theta_zzz = linspace(0,2*pi,500).';
  xxx       = xctot + 3*Rtot*cos(theta_zzz);
  yyy       = yctot + 3*Rtot*sin(theta_zzz);
  zzz = [reshape(xxx,1,numel(xxx));...
	 reshape(yyy,1,numel(yyy))];

  figure(1)
  plot(zzz(1,:),zzz(2,:),'k')		
end
% Assemble the matrix
matAii = get_A_single_diag(C,1:n,kh);

% Build an outer contour
[D,h_D]          = get_circle(contours,1,kh,ppl,1.);
figure(1)
plot(D(1,:),D(2,:),'r')

[Is,Cs,V,matP,D] = get_skeleton(kh,contours{1},matAii,acc,D,h_D);

% Incoming wave
v     = -incoming_wave(kh,theta_in,Ctot(1,:),Ctot(4,:)); v = v.';

% Direct solution
sigma = matAii\v;
phi_scattered_ref  = evalpot(zzz,Ctot,sigma,kh);

% Approximate solution
sigma = matP*v(Is);
phi_scattered      = evalpot(zzz,Cs,sigma,kh);

err = phi_scattered - phi_scattered_ref;

if verif2d
  step     = 4;
  [Inew,idx] = sort(Is);
  Cplot      = Cs(:,idx);
  Cplot      = Cplot(:,1:step:end);
  
  figure(2), hold all;
  err(ind) = NaN;
  err      = reshape(err,ngrid,ngrid);
  contourf(xxx,yyy,log10(abs(err)),-16:2:0)
  cblabel = colorbar('vert');
  plot(Cplot(1,:),Cplot(4,:),'b.','markersize',20)
  plot(D(1,:),D(2,:),'r','linew',2)
  axis equal
  alabel = get(gcf,'CurrentAxes');
  % tlabel = title('outgoing field');
  set(alabel,'FontName','cmr10','FontSize',15,'xscale','lin','yscale','lin','zscale','lin');
  set(cblabel, 'FontName', 'cmr10', 'FontSize', 15);
  % set(tlabel,'FontName', 'cmr10','FontSize',15)
  print(gcf, '-depsc2', '1obj_err.eps');

  
  phi_in		= incoming_wave(kh,theta_in,xxx,yyy);
  phi_scattered_ref	= reshape(phi_scattered_ref,ngrid,ngrid);
  phi_tot		= phi_in + phi_scattered_ref;
  phi_tot(ind)		= NaN;
  figure(3), hold all;
  contourf(xxx,yyy,real(phi_tot))
  cblabel = colorbar('vert');
  plot(Cplot(1,:),Cplot(4,:),'r.','markersize',20)
  axis equal
  alabel = get(gcf,'CurrentAxes');
  % tlabel = title('outgoing field');
  set(alabel,'FontName','cmr10','FontSize',15,'xscale','lin','yscale','lin','zscale','lin');
  set(cblabel, 'FontName', 'cmr10', 'FontSize', 15);
  % set(tlabel,'FontName', 'cmr10','FontSize',15)
    print(gcf, '-depsc2', '1obj_phitot.eps');
else
  fprintf(1,'\n\n ERROR :  %g \n',max(abs(err)))
end

return
